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ABSTRACT 

We explore the use of the bispectrum for understanding quasiperiodic oscillations. The 
bispectrum is a statistic which probes the relations between the relative phases of the Fourier 
spectrum at different frequencies. The use of the bispectrum allows us to break the degen- 
eracies between different models for time series which produce identical power spectra. We 
look at data from several observations of GRS 1915+105 when the source shows strong quasi- 
periodic oscillations and strong broadband noise components in its power spectrum. We show 
that, despite strong similarities in the power spectrum, the bispectra can differ strongly. In all 
cases, there are frequency ranges where the bicoherence, a measure of nonlinearity, is strong 
for frequencies involving the frequency of the quasi-periodic oscillations, indicating that the 
quasi-periodic oscillations are coupled to the noise components, rather than being generated 
independently. We compare the bicoherences from the data to simple models, finding some 
qualitative similarities. 

Key words: accretion, accretion discs - methods: statistical - X-rays:binaries - 
stars:winds,outflows - stars:binaries:close 



1 INTRODUCTION 

Recent studies of the variability properties of accreting black holes 
and neutron stars have shown that these systems' Fourier power 
spectra in certain spectral states are well-described by a sum of 
many Lorentzian components (see e.g. Olive et al. 1998; Nowak 
2000; Belloni, Psaltis & van der Klis 2002; van Straaten et al. 2002; 
Pottschmidt et al. 2003). Essentially the same components seem to 
appear in nearly all sources, and at nearly all luminosities, and their 
frequencies tend to be well correlated with one another and with the 
source spectral states (Wijnands & van der Klis 1999; Psaltis, Bel- 
loni & van der Klis 1999). This phenomenology suggests that there 
may be a single origin for most of the variability features in the 
power spectra of accreting black holes and neutron stars. Similar 
correlations in frequencies between the different components seem 
to apply even to accreting white dwarfs (e.g. Mauche 2002; Warner, 
Woudt & Pretorius 2003). 

Coupling between different components in the power spec- 
trum has been suggested in many theoretical contexts. In a shot 
noise model (e.g. Terrell 1972), variability components on the 
timescale of the shot will be correlated with one another, result- 
ing in, e.g. fast rise exponential decay or exponential rise, fast de- 
cay profiles, but there should be no coupling on timescales longer 



than those of the shots. More sophisticated models of variability, 
for example, self-organized criticality (e.g. Takeuchi, Mineshige & 
Negoro 1995; Takeuchi & Mineshige 1997), predict correlations 
between the arrival times and/or intensities of the shots. Other re- 
lated "reservoir" models (e.g. Merloni & Fabian 2001; Maccarone 
& Coppi 2002; Malzac, Merloni & Fabian 2004), where the ac- 
cretion disk and/or the relativistic jet taps an energy supply effec- 
tively enough to reduce the available energy for future emission 
also predict variability correlated over many frequencies. Reso- 
nance models for producing quasi-periodic oscillations should also 
clearly produce non-linear variability (e.g. Psaltis & Norman 2001; 
Abramowicz & Kluzniak 2001; Schnittman & Bertschinger 2003; 
Maccarone & Schnittman 2005), while non-resonant models for 
producing the same QPOs (e.g. Chen & Taam 1992,1995; Rezzolla 
et al. 2003; Giannios & Spruit 2004) could, but need not show cou- 
pling between the different frequencies. Propagation models, where 
disturbances move through an accretion disc, also should produce 
non-linear coupling of variability components, as the variability 
properties are modified at each annulus (Lyubarskii 1997). 

A key to verifying and understanding this possible unified ori- 
gin for variability is to go beyond the simple power spectrum and 
begin studying their non-linear variability. The first attempts at this 
failed to detect signatures of non-linearity, with the methods includ- 
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ing the time skewness (Priedhorsky et al. 1978) and searches for a 
low dimensional chaotic attractor (Lochner, Swank & Szymkowiak 
1989). On the other hand, more recent work with better light curves 
did establish that the light curves of Cygnus X-1 are not time re- 
versible (Timmer et al. 2000; Maccarone & Coppi 2002), that there 
exists a correlation between rms amplitude and flux of a source that 
is inconsistent with pure shot noise models (Uttley & McHardy 
2001), that there is coupling between variability components on 
all observable time scales (Maccarone & Coppi 2002), and that in 
some cases, there is a low dimensional chaotic attractor after all 
(Misra et al. 2004), and the observed light curves may be described 
by a Lorenz system (Misra et al. 2006). With this in mind, we now 
approach looking at the properties of the coupling between quasi- 
periodic components and noise components in GRS 1915+105, a 
bright Galactic X-ray binary with an accreting black hole which 
has been the subject of several long observations with the Rossi 
X-ray Timing Explorer (RXTE). 

We treat this work as a pilot study - the first attempt to ap- 
ply the bispectrum to astronomical data with strong quasi-periodic 
oscillations. As such, we aim to illustrate the power of the tech- 
nique by showing data and simple models that can give very similar 
power spectra, and very different bicoherences, but consider it be- 
yond the scope of the work to try to fit models to the data precisely. 
We present computations of the bicoherence for several observa- 
tions of GRS 1915+105, discuss the meaning of the bicoherence, 
and present a few toy models for the bicoherence which we com- 
pare with the data. 



2 DATA USED 

For this pilot study, we present the results from three representa- 
tive observations of GRS 1915+105; our goals in this paper are not 
to make a complete characterisation of the bicoherence properties 
of all sources in all states, but rather to demonstrate of the util- 
ity of the technique for identifying phenomenological differences 
between apparently rather similar lightcurves with rather similar 
power spectra, and to show a few theoretical models that produce 
bicoherences similar to some of the observations. All these data are 
from variability class x in the classification scheme of Belloni et 
al. (2000). In this class, GRS 1915+105 stays steadily in state C, 
its "low/hard" state - the "hard very high state" in the nomencla- 
ture of Fender & Belloni (2004), or the hard intermediate state in 
the more recent classification of Homan & Belloni (2005). In this 
state, GRS 1915+105 shows a strong power law component in its 
X-ray spectrum, and high rms amplitude quasi-periodic oscillations 
at about 0.5-10 Hz, often with powerful harmonics (see Fender & 
Belloni 2004 and references within). The broadband noise compo- 
nent in the power spectrum containts a greater fraction of the rms 
amplitude of variability than do the QPOs and harmonics. 

The RXTE observation identification numbers of the three 
observations presented are 10408-01-25-00 (taken 19 July 1996), 
20402-01-15-00 (taken 9 Feb 1997), and 30184-01-01-000 (taken 
4 April 1998). Standard screening has been applied to the data to 
ensure no usage of data taken during Earth occulations or periods 
of high offset. The data used are the single bit mode data with the 
lowest set of energy channels: channels 0-35 are used for 10408- 
01-25-00 and 30184-01-01-000, while channels 0-13 are used for 
20402-01-15-00. These correspond approximately to the energy 
range from 2-13 keV and 2-5 keV, respectively. We have done 
checks of the other data modes with different energy channel ranges 
and have found no qualitative differences and only minor quantita- 



tive differences in the properties of the bicoherence as a function 
of photon energy, so we have used the lowest energy set of chan- 
nels in a single data mode for each observation for convenience - 
we will not attempt to make interpretations of the data at a level of 
detail approaching the differences caused by changing energy band 
used. The data are then Fourier transformed with 4096 elements per 
transform, and time resolution of 1/64 seconds, 1/128 and 1/512 
seconds, respectively. These are chosen so that the quasi-periodic 
oscillation is near the middle of the range of Fourier bins in each 
observation (i.e. a lower time resolution is used when the QPO has a 
lower frequency). The integration times used are 8640, 10240, and 
8160 seconds, respectively, with only the first part of observation 
30184-01-01-000 used. This decision was originally made due to 
memory limitations when the calculations were first made, and the 
decision not to re-do the calculation was made because there is sig- 
nificant frequency drift over the full observation, and the effects of 
frequency drift are ameliorated by using shorter total integrations. 
The Fourier transforms are then combined to produce the bispectra 
and bicoherences of these data sets. 

We have examined the amount of frequency drift in the dif- 
ferent observations. For observation 10408-01-25-00 and 20402- 
01-15-00, the peak frequency in the power spectrum varies by less 
than 10% over the full integrations - less than the widths of the 
QPOs. For 30184-01-01-000, there is substantial variation of the 
peak frequency if we consider the full observation - with peak fre- 
quencies ranging from about 2.7 Hz to 4 Hz. However, in the part 
of the data set we consider, the variation is only from about 3.3 Hz 
to 4 Hz. Some minor effects from this variation can be seen in the 
computed bicoherences, and we discuss them below. 



3 STATISTICAL METHODS 

For this work, we will focus on computing the bispectrum and the 
closely related bicoherence. The bispectrum computed from a time 
series broken into K segments of equal length is defined as: 

K-l 

B(/c,0 = ^Xi(fc)Xi(OX*(/c-FZ), (1) 

i=0 

where Xi{k) is the frequency k component of the discrete Fourier 
transform of the ith time series (e.g. Mendel 1991; Fackrell 1996 
and references within). The bispectrum is well defined only for 
k + I less than or equal to the Nyquist frequency of the data used 
to compute it, and is defined nontrivially only for k < I, since 
B{k, I) = B{1, k) - although in this paper, we will plot the bico- 
herence for fc > / and for / > fe to help guide the eye to features 
in the plots. The expectation value of the bispectrum is unaffected 
by additive Gaussian noise, although its variance will increase for a 
noisy signal. Poisson noise can affect the bicoherence for frequen- 
cies where the Poisson level is high compared with the level of 
intrinsic variability - the amplitude of the Poisson noise correlates 
with the count rate, so that there will be phase coupling between the 
instrinsic variability of the source's signal, and the noise level in the 
Poisson component (see e.g. Appendix C of Uttley et al. 2005). In 
this paper, we focus on low frequency variabiUty of a highly vari- 
able source, so Poisson noise is unimportant. 

The definition of the bispectrum gives its absolute phase a 
well-defined meaning, in constrast to the phase of the Fourier spec- 
trum which has a dependence on an arbitary reference time. One 
can then attempt to determine the degree of constancy of that phase, 
making use of a related quantity, the bicoherence - the magnitude 
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of the bispectrum, normalised to lie between and 1 . Defined anal- 
ogously to the cross-coherence function (e.g. Nowak & Vaughan 
1996), it is the vector sum of a series of bispectrum measurements 
divided by the sum of the magnitudes of the individual measure- 
ments. If the biphase (the phase of the bispectrum) remains constant 
over time, then the bicoherence will have a value of unity, while if 
the phase is random, then the bicoherence will approach zero in 
the limit of an infinite number of measurements. Specifically, the 
squared bicoherence, is defined as: 



\Y,Mk)X^[l)X*{k + l)\ 

Y.\x,{k)x,{i)\'Y.\Mk + i)? 



(2) 



This normalization of the bicoherence was proposed by Kim 
& Powers (1979). If it shows variation as a function of the frequen- 
cies used, the time series analysed is nonlinear. It has been shown 
that other normalisations are more likely to detect nonlinearity un- 
der certain specific circumstances (Hinich & Wolinsky 2004), but 
we use the Kim & Powers (1979) normalisation here to retain con- 
sistency with past work. In previous work where we were studying 
variability components which were very broad in the power spec- 
trum (Maccarone & Coppi 2002), we binned the data to make a 
one dimensional function of A; + and substantially re-binned the 
bispectrum values over frequencies and then compared the results 
with model predictions. Since that work, we have become aware of 
a correction which is, in principle, important for studying aperiodic 
variability with the bicoherence. The maximum value of the bico- 
herence is suppressed by smearing of many frequencies into a sin- 
gle bin in the discrete Fourier transform (S. Vaughan, private com- 
munication). This suppression cannot be calculated in a straightfor- 
ward way (see e.g. Greb & Rusbridge 1988). We note also, though, 
that since comparisons in Maccarone & Coppi (2002) were made 
only with model calculations made with the same time binning as 
the real data, these effects, whatever they may be, are the same for 
the real data and the simulated data, and hence the conclusions of 
that paper are not affected substantially. Therefore, to make quan- 
titative comparisons of real data with simulated bicoherences, it is 
important to use the same frequency binning in both cases. 

A non-zero bicoherence can be used to rule out, for example, 
models where the variability on different timescales comes from in- 
dependent Gaussian model components which are then added lin- 
early. Beyond that, it is often difficult to interpret the results of 
bicoherence analysis, but comparison of observed results with re- 
sults from simulated light curves can be a very powerful tool. As 
shown already, for example, in Maccarone & Schnittman (2005), 
light curves with very similar power spectra can have easily distin- 
guishable bicoherence plots. 

For this work, we will consider variability only at frequencies 
which are low enough that the intrinsic variability of the source 
is strong, so that neither Poisson noise nor dead time affects the 
Fourier spectrum substantially. We do note that these problems may 
be more important at higher frequencies and that they will need to 
be considered in future analyses. 



4 BICOHERENCE RESULTS 

We have computed the bicoherence for many observations of 
GRS 1915+105. For the sake of brevity, we present only three rep- 
resentative results, aiming to present three observations with very 
similar power spectra, but qualitatively different bicoherences. We 
find three phenomenological patterns of variability here, which we 
described as "web", "cross", and "hypotenuse". The bicoherences 



from actual data are plotted in figure [T] and the corresponding 
power spectra are plotted in figure[2l 



4.1 The "cross" pattern 

When the source is in the "radio quiet" state C at low count rates, 
the QPO frequency is low (less than 2 Hz). Here, the bicoherence 
shows a "cross" behavior (see the results from ObsID 20402-01-15- 
00 plotted in figureflj) for the bicoherence and figure|2j5 for the cor- 
responding power spectrum), with large bicoherence for frequency 
pairs where one frequency is the QPO frequency and the other can 
take on any value. Here, in contrast to the web pattern, there is 
large bicoherence even for noise frequencies less than twice fgpo- 
Additional power is seen at the harmonics and for the harmonics 
interacting with the noise. 

It is interesting to consider whether the cross pattern might 
be indicating the effects of very broad wings to the QPO, with 
no actual coupling between the QPO and the noise, but just cou- 
pling between the QPO and its harmonics. This is unlikely to be the 
case. Firstly, the bicoherence features are asymmetric in frequency 
space, while the QPOs are well-fit by Lorentzians, indicating that 
they are not asymmetric in frequency space. In principle, since the 
bicoherence is normalised by dividing by quantities related to the 
power spectrum, changes in the strength of the noise component, 
which would dilute the coupled variability, could be expected. In 
this case, though, there is both more power in the noise component 
of the power spectrum and stronger bicoherence at frequencies be- 
low the QPO frequency than above it. Secondly, the harmonics tend 
to follow either roughly circular patterns for their bicoherence, if 
the quasi-periodicity is caused by phase disconnections in an oth- 
erwise periodic signal, or /i = /2 patterns, if the quasi-periodicity 
is caused by the presence of many real frequencies due to changes 
in frequency over time or the superposition of multiple frequencies 
as might happen if the variability is due to orbital motions over a 
range of radii (see Maccarone & Schnittman 2005), and what is 
observed is different from both. 



4.2 The "hypotenuse" pattern 

At higher count rates of the "radio quiet" state C, the bicoherence 
pattern changes dramatically (see the results from ObsID 30184- 
01-01-000 plotted in figure[T}; for the bicoherence and figure|2j; for 
the corresponding power spectrum) . A high bicoherence is seen 
primarily where the two noise frequencies add up to the QPO fre- 
quency - that is, the regions of large bicoherence make a diagonal 
line. We refer to this as the "hypotenuse" pattern because the re- 
gion of strong bicoherence forms the hypotenuse of a triangle with 
its other two sides being the axes of the plot. The "hypotenuse" 
pattern typically also shows a roughly circular spot of high bicoher- 
ence with /i=/2=/qpo, since there is power both at the fundamen- 
tal and at the first harmonic and these two variability components 
are coupled. A few other observations at similar count rates and 
with similar QPO frequencies in the non-plateau x state have been 
examined, and show the same qualitative features in their bicoher- 
ences. We also note that the diagonal elongation in the bicoherence 
plot from lower left to upper right for fi = f2 ~ Jqpo is likely 
to be related to the frequency drift during the observation that is 
mentioned above. 
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4.3 The "web" pattern 

A hybrid class of bicoherence plots can be seen for some of the ob- 
servations, where features resembling those seen in the "cross" and 
of the "hypotenuse" patterns can be seen together. When the source 
is in the radio plateau state, as in ObsID 10408-01-25-00, the bi- 
coherence shows a rather strong signal for two combinations (see 
figure [TJi and figure for the corresponding power spectrum) - a 
diagonal line from upper left to lower right, with /i + /2 = /qpo, 
for all < /i,/2 < /qpo- Additionally, one vertical/horizontal 
streak is seen with one of the frequencies equal to the QPO fre- 
quency, and the other a frequency greater than the twice QPO fre- 
quency. No significant vertical/horizontal signal is seen for the case 
that /i < 2/qpo. The value of the bicoherence tapers off for fre- 
quencies where the power in the noise component becomes small. 
We refer to this pattern of behavior as the "web" pattern. A few 
other observations at similar count rates in the plateau state have 
been examined, and show the same qualitative features in their bi- 
coherences. 



4.4 A few brief remarks on other observations 

In some observations not plotted in this paper, with the highest 
count rates seen from GRS 1915+105 in its x class, the bicoher- 
ence seems to tend towards zero entirely. However, at the highest 
count rates, the QPO frequency can also change substantially on 
the timescale of a single RXTE observation. The bicoherence is 
not well suited to dealing with data sets which are non-stationary 
in this manner. This problem may exist even in some of the lower 
count rate data sets, but at a much lower level. We note that in a 
few of the observations, the bicoherence is extended along the line 
fi ~ f2 around the location where the interactions are between 
the QPO and its first harmonic - this effect can be seen showing 
up weakly in Figure la. It was shown in Maccarone & Schnittman 
(2005) that this is characteristic of a QPO which is broadened due 
to the existence of power at many frequencies, rather than due to 
phase disconnections. 



5 DISCUSSION 

The overall values of the bicoherence in these data can also give 
us some clues as to the types of processes that might be producing 
the phase correlations. The bicoherence relates the instantaneous 
response of the driven oscillation to changes in the driving modes. 
As stated above, Greb & Rusbridge (1988) have shown that the 
maximum value of the squared bicoherence is Suj/T,, wheieSuj is 
the frequency resolution in the observation (or the width of a QPO, 
if the QPO is resolved), and E is the frequency width of the driving 
spectnmi. Therefore, if a QPO is driven by interactions of all noise 
frequencies less than the QPO frequency, then the largest possible 
value of the squared bicoherence is Q~^, where Q is the quality 
factor of the QPO. 

Greb & Rusbridge (1988) also consider the effects of having a 
relatively narrow resonance. In this case, the squared bicoherence is 
reduced by a factor of the ratio of the coherence time in the driven 
mode to that in the driving mode. Thus, for resonant modes, which 
remain coherent for long periods of time relative to the timescale 
on which the driving signal remains coherent, the bicoherence will 
be substantially reduced. This can be grasped intuitively in a quali- 
tative sense - a strong resonance will have an amplitude related to 



the integral of the power dumped into it over a long timescale, and 
hence not strongly correlated with the instantaneous driving power. 

In our observations, the Q values of the QPOs are typically 
~ 3—10, and the QPOs are overresolved in frequency space by fac- 
tors of about 10-20. This yields expected maxima for the squared 
bicoherence of ~ 0.05 (with the exact value depending both on 
which observation is considered, and whether multiple noise com- 
ponents interact nonlinearly to produce the QPO, or the QPO in- 
teracts nonlinearly with one noise component to produce another 
noise component). Larger values are possible for the harmonics of 
the QPOs, where the width of the driving spectrum is smaller than 
the width of the noise component in frequency space. 

If we assume that the perturbations which are coupled to one 
another give an X-ray count rate which is linearly proportional to 
the amplitude of the perturbation, we can then use the magnitude of 
the bicoherence to gain some insight into whether the interactions 
can be through a narrow resonance at the QPO frequency. Our find- 
ings that the maximum values of the squared bicoherence for the 
cases where the QPO and the noise are interacting are ~ 10~^'^ 
would then imply that the coupling between the noise and QPO 
is relatively close to its maximal value. This places some immedi- 
ate constraints on the classes of models which can be considered 
- models must have a large fraction of the power coming from an 
emission region which behaves as a single "system", and, if the 
interactions between noise components produce a QPO through a 
resonant interaction, it must be a highly damped resonance. 

5.1 Simulated light curve analysis 

We now consider several different mathematical forms for light 
curves which correspond, at least approximately, to physical sce- 
narios for producing QPOs and noise components which are cor- 
related in some manner. We show that these different models can 
produce relatively similar power spectra while producing bicoher- 
ence diagrams which are qualitatively quite different from one an- 
other. Two of the patterns of behaviour found - the "cross" and 
the "hypotenuse" - are reproduced reasonably well by relatively 
straightforward models. The third, the "web" is only approximately 
reproduced here with a relatively simple model. 

5.1.1 Bicoherence in terms of reservoir models 

Where the bicoherence's value is large, the variability must be cou- 
pled on the three timescales corresponding to the three frequen- 
cies included for the computation of the bicoherence. A natural 
way to couple variability on different timescales is with a reservoir 
model, where there are multiple variability components that drain 
the same reservoir. In this way, the different components "compete" 
for power - when the noise component becomes stronger, there will 
be less energy in the reservoir available for the QPO and vice versa, 
leading to a phase coupling. An analogous idea has been put forth 
to explain the coupling between optical variability and X-ray vari- 
ability in XTE Jl 1 18+480 (Malzac, Merloni & Fabian 2004). 

To determine what the output light curve will look like for the 
simulations, we define, in each case, two power spectra to use as 
outputs. One, the noise component, is a broad Lorentzian with a 
peak frequency of zero. The other, the quasi-periodic oscillation, is 
a much narrower Lorentzian with a peak at a higher frequency (we 
use Q = 30 here, where Q is the quality factor of the Lorentzian). 
The time series for these components are made using the method 
of Timmer & Konig (1995). This method generates simulated times 
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series from a power spectrum under the assumption that the process 
is Gaussian - thus we can be assured that any non-linearity we de- 
tect is because of the physics we put in after generating the initial 
random time series. 

Then, a reservoir is defined. Energy is injected into the reser- 
voir either at a constant rate, or at a purely random rate (i.e. a ran- 
dom number uniformly distributed between zero and one times a 
normalisation). Energy is drawn from the reservoir by each of the 
two components. The flux of a component is taken to be the value 
of its time series multiplied by the size of the reservoir at that in- 
stant times some normalization. The flux is then subtracted from 
the reservoir. Specifically, we do the following: 

R'(t)=R{t-M)+Yr{t), (3) 

where R' it) is the size of the reservoir at time t, after the energy 
injection has taken place, but before the draining of the reservoir 
has taken place, At is the time step, R(t — At) is the size of the 
energy reservoir at the end of the previous time step, and Yr{t) is 
the injection rate into the reservoir per time step; 

Ygpoit) = ygpoit) x R'{t) x kQPo, (4) 

where Vqpo is the output flux from the QPO component, yqpo 
is the time series of the QPO produced by the Timmer & Konig 
method, and kqpo is a normalization constant; 

YN{t) = yN{t) X 7?'(t) X fcjv, (5) 

where Yn is the output flux from the noise component, i/jv is the 
time series of the noise produced by the Timmer & Konig method, 
and fcjv is a normalization constant; and, finally, 

R{t) = R'{t) - YQPo{t - At) -YN{t- At). (6) 

Two general cases are considered, and we search a range of 
parameter space for each case. The first is where the output "light 
curve" is the sum of the fluxes of the two components. This might 
correspond to a case where accretion power is dissipated either in 
some oscillating region which produces the QPO, or a non-resonant 
part of the accretion disc, which produces the noise component. 
The second is where the output light curve is only the QPO's flux. 
This might correspond to the case where the broad Lorentzian rep- 
resents a jet which extracts power from the system, but does not 
emit in the X-rays (see e.g. Malzac et al. 2004). 

If the reservoir is drained by a component which does not con- 
tribute to the X-ray light curve, then the size of the reservoir in 
"view" of the X-ray light curve, will effectively be some quantity 
minus the integral of the power drained by the other component. 
As a result, the reservoir itself will have a signature of its modula- 
tion, so when the output light curve is obtained by multiplying the 
reservoir size by the time series for the QPO, the "QPO" compo- 
nent will now be modulated on the timescale of the "noise," and the 
noise power spectrum will affect final power spectrum. As we show 
below, the qualitative properties of the bicoherence are largely the 
same regardless of whether the two drains on the energy reservoir 
both contribute to the output flux, or only one of these does. Again, 
the logic we follow here is quite similar to that in Malzac et al. 
(2004). Given that, as we will show below, these models seem to 
match the data for the plateau states more closely than the more 
radio quiet x states, it would not be surprising if energy extraction 
by a jet were an important factor in determining the properties of 
the light curve. 

After the basic model equations are defined, a time series is 
produced for the model. The QPO and noise components are pro- 
duced as time series using the Timmer & Konig method, with both 



components assumed to be Lorentzians with parameter values set 
to values similar to those seen in the data for various observations. 
For a very small number of time bins, the value of one of these 
time series will be negative, in which case we set it to zero. Alter- 
natively, the energy extracted in a time interval may be larger than 
the size of the reservoir, in which case we set the power such that 
the whole reservoir is drained in that time step. 

When the amplitude of variability is large, the reservoir's level 
fluctuates strongly. This leads to the production of harmonics in the 
QPO (since the oscillations become non-sinusoidal when the reser- 
voir fluctuates in response to the oscillations themselves). It also 
leads to phase coupling of all types between the QPO and the noise. 
This produces a "web-like" pattern in the bicoherence; however, the 
web-like pattern here is not identical to the one seen in the radio 
plateau state. In the real data, the "cross-like" structures begin to 
manifest themselves only for noise frequencies larger than the fre- 
quency of the first harmonic, while in these simulations, the cross 
structures appear at all noise frequencies. Furthermore, in some of 
the simulations, there are diagonals from upper left to lower right 
related to the cross structure which are strong features at all fre- 
quencies adding up to at least one of the harmonics. These are 
not present in the real data. The strong diagonals manifest them- 
selves when the total variability in the reservoir is dominated by 
the noise component, rather than by the QPO. The model bicoher- 
ence is shown in Figure 3b to represent a case where the jet extracts 
power (i.e. so that the extracted noise-component power does not 
lead directly to observed emission), and in Figure 4b for the case 
where the noise component is added back in to the model light 
curve. Figure 5b shows the simulated power spectrum correspond- 
ing to the bicoherence shown in Figure 3b. 

In the numerical calculations where the noise component is 
not added to the final light curve, the QPO has a mean count rate 
of 8000, with a root mean squared variability level of 1200, while 
the noise component has a mean value of 20000, with an root mean 
squared variability level of 8000. We re-fill the reservoir by adding 
a random number drawn from a uniform distribution between and 
8000. The draining of the reservoir takes place by summing the 
two time series, multiplying by the size of the reservoir, dividing 
200000, and subtracting that from the reservoir value. 

Where the noise component is added to the final light curve, a 
slightly different procedure is adopted. The QPO and noise compo- 
nents are added together in the initial power spectrum made before 
the inversion into a time series using the Timmer & Konig method. 
The QPO is given a normalization in the power spectrum 200 times 
higher than the noise component's. The actual value in the final 
power spectrum comes from converting this power spectrum into 
a time series with mean value of 8000 and rms amplitude of 6000. 
All other procedures are as above. 

At low variability amplitude, only the cross-like structures are 
seen strongly, although a weak "hypotenuse" feature is still present 
when the noise component is added to the output "light curve." The 
bicoherence is shown without the noise component added back in 
in Figure 3a, and with it added back in in Figure 4a. The power 
spectrum corresponding to the low variability amplitude model is 
shown in Figure 5a. The simulated bicoherence plots in this case 
look quite a bit like the bicoherence plot seen from observation 
20402-01-15-00 (compare the observational data in FigureflJ) with 
the simulations in Figure[3^). 

The procedure for converting the model into a simulated light 
curve is the same as above, except for some changes in values of 
parameters. For the case where the noise component is not seen 
in the output simulated light curve, the noise component in this 
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case has an rms amplitude of 1000 instead of 8000, with all other 
parameter values the same. Where the noise component is added to 
the final simulated light curve, the rms amplitude of the simulated 
light curve is 2500 instead of 6000. 

5.1.2 Bicoherence in terms of damped forced harmonic 
oscillators 

We have also considered the case of a damped, forced harmonic os- 
cillator. This is conceptually similar to the mechanism suggested by 
Psaltis &. Norman (2000) for producing both quasi-periodic oscil- 
lations and noise components from a single perturbation spectrum 
- the quasi-periodic oscillation appears at the resonant frequency 
of the system, while a white noise input spectrum is turned into a 
red noise component by the fact that damped oscillators act as low- 
pass filters. However, the response of the system was considered to 
be linear in that work, meaning that non-linear features like non- 
zero bicoherence can not be explained without some modifications. 
When modifying the system of equations, to allow for a non-linear 
restoring force, we find that the resulting pattern of bicoherence that 
results from this scenario is relatively similar to the "hypotenuse" 
pattern, although we do have trouble reproducing a strong bico- 
herence with a realistic power spectrum. An alternative physical 
mechanism for producing such a quasi-periodic oscillation would 
be to have variations in the accretion rate excite inertial-acoustic 
modes in the inner accretion disc (Chen & Taam 1992; 1995). 

Quantitatively, this scenario is modelled in the same way as a 
spring which does not follow Hooke's Law (i.e. F — ~kx), but 
rather has an asymmetric, higher order restoring force, F which is 
given hy F = ~kix + k2x'^, where x is the displacement, and fci 
and k2 are constants parameterizing the strength of the restoring 
force. If k2 — 0, this does reduce to Hooke's Law, and the lack 
of a non-linear term leads to a lack of observable features in the 
bicoherence plot. 

The full equations of motion for this "oscillator" are: 

F{t) = -kix{t) + k2X^{t) - "fv{t) + YwN{t), (7) 

where 7 is the damping coefficient and Ywn is the magnitude of 
the driving force, which is assumed to be a white noise process; 

v{t) = v(t- At) + At X F(t) /m, (8) 

where m is the mass of the oscillator, which we set to unity for these 
calculations, and At is 0.001. The calculated simulated light curve 
is the time series x{t), and has 1048576 points. The bicoherence is 
computed from 256 chunks of 4096 points. 

For the first calculation, we set the parameter values to ki = 
50000.0, k2 = 0.2, and 7 = 7.0. We show the plot of simulated 
data in[3]:. We note there is one substantial difference between the 
model data and the real data for 30184-01-01-000, which is that the 
model data has stronger bicoherence along the "hypotenuse" when 
the two frequencies are most different from one another, while the 
real data shows the opposite trend. This trend makes the model 
values more seriously deviant from the real data for the strongest 
damping of the motion of the oscillator. 

When the damping factor is very low, then the resulting bico- 
herence plots can look quite similar to the "web" pattern. In figure 
[6] we show the power spectrum and bicoherence plot for a low 
damping case. The calculations are done in the same manner as the 
ones described in the previous paragraph, except with ^2=0.25 and 
7 = 0.6. For values of 7 intermediate between 0.6 and 7.0, the "hy- 
potenuse" can become much broader without creating such strong 
additional harmonics in the power spectrum. This simulation, like 



the one whose power spectrum and bicoherence are shown in Fig- 
ures 4213 and[3j3, respectively, matches many of the properties of 
the data, but not all; however, here, the problem is primarily that 
the model shows too strong a bicoherence in streaks where there is 
only noise, rather than power in the QPO or its harmonics. 

5.1.3 Caveats 

There are further steps which will be needed to make models which 
match all the statistical properties of the observed light curves. For 
example, in most cases, the observed flux distributions will follow a 
log-normal distribution (Uttley et al. 2005) - we do not attempt here 
to enforce this requirement. We do note that it was found by Utt- 
ley et al. (2005) that exponentiating the flux values in a light curve 
can push a flux distribution towards a log-normal distribution. We 
have done this in a few cases, and have found that the bicoherence 
plots change very little when the flux values are exponentiated. Ut- 
tley et al. (2005) have already shown that power spectra are not 
typically affected strongly by this transformation. Therefore, flux 
distributions and bicoherence plots are largely independent tests of 
the variability's coupling, although the flux distributions produced 
in these simulated data sets do not deviate strongly from log-normal 
distributions in any event. Furthermore, in cases where there are 
strong QPOs with very low frequencies, there are clear deviations 
of the flux distribution from being log-normal (Misra et al. 2006). 
We have additionally not ensured that the power spectra of real ob- 
servations are matched exactly by these simulations. In essence, 
since the input models do not have any radiative transfer in them, 
we are assuming that the radiative efficiency of any perturbation 
will be constant when we call the output time series "simulated 
light curves." Deviations from this assumption are likely to affect 
the power spectra, and, especially, the flux distributions, far more 
seriously than they affect the morphology of the bicoherence plots. 
Such deviations are quite easy to imagine in situations where there 
exists a resonance. 

An additional potential problem exists with connecting the toy 
models presented here to reality. Much of the phenomenology ap- 
pears to be similar between the different sub-types of x classes. It 
would therefore be surprising if there were a qualitative change in 
the mechanism for producing the QPOs. It would thus be useful to 
identify a means to unite the different mechanisms described here, 
or to develop an entirely new physical model to explain all the ob- 
served phenomenology in a more unified way, with just the change 
of one or two parameters causing the broad range of phenomena 
seen in terms of bicoherence patterns. 



6 CONCLUSIONS 

We have shown clearly that the variability in the noise compo- 
nents and quasi-periodic oscillations of GRS 1915+105 are cor- 
related with one another. The bicoherence provides a good dis- 
criminator among variability models which produce similar power 
spectra through quite different physical processes. We have found 
that it is plausible that in the "plateau" state of GRS I9I5+I05, 
the variability is caused by a reservoir of energy being drained by a 
noise component (which could be the radio jet) and a quasi-periodic 
component, while in the brighter part of the x state, the variability 
is consistent with a white noise input spectrum driving a damped 
harmonic oscillator with a non-linear restoring force. While the 
models presented here are almost certainly not unique solutions for 
what is occuring physically in these systems, it is quite clear that 
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Figure 1. The bicoherence plots from the observations, (a) the data for 
observation 10408-01-25-00, the prototypical "web" source. The colour 
scheme is given for different values of the log^gb^ as follows: red: -1.0; 
green: -1.25; dark blue: -1.5; light blue: -1.75 (b) the data for observation 
20402-01-15-00, the prototyical "cross" source. The colour scheme is given 
for different values of the logj^Qft^ as follows: red: -1.0; green: -1.25; dark 
blue: -1.5; Hght blue: -1.75 (c) the data for observation 30184-01-01-000, 
the prototypical "hypotenuse" source. The colour scheme is given for dif- 
ferent values of the logj^gfe'^ as follows: red: -1.0; green: -1.25; dark blue: 
-1.5; light blue: -1.75. The equality of the values of the bicoherence for 
reflections about x = yif, trivial. 
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Figure 2. The power spectra for the real data: (a) observation 10408-01-25- 
00 (b) observation 20402-01-15-00 (c) observation 30184-01-01-000. 



the bicoherence will provide excellent constraints on more physi- 
cally motivated models for the variability, and that we can definitely 
identify cases where the properties of the power spectra are quite 
similar, but the properties of the light curves are quite different. 
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Figure 3. The model bicoherence plots. For all plots, the colour schemes 
are as for the real data, except that for the model plots, there is an addi- 
tional contour, in purple, at logj^pfe^ = —2.0. The plots are: (a) the case 
of "jet" extraction with a high variability amplitude (b) the case of "jef 
extraction with a low variability amplitude (c) damped forced oscillations. 
The equality of the values of the bicoherence for reflections about x = yis 
trivial. 
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Figure 4. The model bicoherence plots for reservoir models where the 
noise component contributes to the observed flux. For all plots, the colour 
schemes are as for the previous plot of simulated data. The plots are: (a) 
low variability amplitude (b) high variability amplitude. The equality of the 
values of the bicoherence for reflections about x = y is trivial. 
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Figure 5. The power spectra for the simulated data: (a) web-hke (b) cross- 
Uke (c) damped forced oscillations - note that these are the power spectra 
for the same time series, in the same order, as presented in Figure 3. 



